WEAK GALERKIN METHODS FOR SECOND ORDER ELLIPTIC 
INTERFACE PROBLEMS 
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Abstract. Weak Galerkin methods refer to general finite element methods for partial differ- 
ential equations (PDEs) in which differential operators are approximated by their weak forms as 
distributions. Such weak forms give rise to desirable flexibilities in enforcing boundary and inter- 
face conditions. A weak Galerkin finite element method (WG-FEM) is developed in this paper for 
solving elliptic PDEs with discontinuous coefficients and interfaces. The paper also presents many 
numerical tests for validating the WG-FEM for solving second order elliptic interface problems. For 
such interface problems, the solution possesses a certain singularity due to the nonsmoothness of 
C^^ the interface. A challenge in research is to design high order numerical methods that work well for 

fNj problems with low regularity in the solution. The best known numerical scheme in the literature is 

of order 0(h) for the solution itself in Loo norm. It is demonstrated that the WG-FEM of lowest 
order is capable of delivering numerical approximations that are of order 0(h^-^^) in the usual Loo 
norm for C-^ or Lipschitz continuous interfaces associated with a or continuous solution. The- 
oretically, it is proved that high order of numerical schemes can be designed by using the WG-FEM 
with polynomials of high order on each element. 
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1. Introduction. In this paper, we are concerned with interface problems for 
second order elhptic partial differential equations (PDEs) with discontinuous coeffi- 
cients and singular sources. For simplicity, consider the model problem of seeking 
functions u = u{x^y) and v = v{x,y) satisfying 



(1.1) 


- V.AiVi/ = /i, 


in r^i. 


(1.2) 


-V'A2Vv = f2, 


in 1^2, 


(1.3) 


u = gi, 


on dQi \ r. 


(1.4) 


V = ^2, 


on \ r. 


(1.5) 


u — V = (j), 


on r. 


(1.6) 


AlVu • ni + A2\/v • n2 = '0, 


on r. 
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.J^ normals of and Q2- Here, /i and /2 can be singular and F may be of Lipschitz 
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continuous. This problem is commonly referred to as an elliptic interface problem and 
occurs widely in practical applications, such as fluid mechanics [23|, electromagnetic 
wave propagation [Mj [181 Hll HO] , materials science ^Ol |22] , and biological science 
[39l [T2l The finite difference based solution of elliptic interface problems was 
pioneered by Peskin with his immersed boundary method (IBM) in 1977 [30l |29] . 
Mayo constructed an interesting integral equation approach to this class of problems 
[25 . To properly solve the elliptic interface problem, one needs to enforce additional 
interface conditions (1.5) and (1.6). LeVeque and Li advanced the subject with their 
second order sharp interface scheme, the immersed interface method (JIM) [24 . In the 
past decades, many other elegant methods have been proposed, including the ghost 
fluid method (GFM) proposed by Fedkiw, Osher and coworkers [11], flnite- volume- 
based methods [28 , the piecewise-polynomial discretization [HI , and matched interface 
and boundary (MIB) method ^[43l|37]. 

A proof of second order convergence of the IIM for smooth interfaces was due to 
Beale and Layton [2 . Rigorous convergence analysis of most other flnite difference 
based elliptic interface schemes is not available yet. In general, it is quite difficult 
to analyze the convergence of flnite difference based interface schemes because con- 
ventional techniques used in Galerkin formulations are not applicable for collocation 
schemes. The analysis becomes particularly difficult when the designed elliptic inter- 
face scheme is capable of dealing with nonsmooth interfaces [37]. At present, there is 
no rigorous convergence analysis available for elliptic interface methods that deliveries 
high-order accuracy for nonsmooth interfaces, to the author's knowledge. 

Finite element methods (FEMs) are another class of important approaches for 
elliptic interface problems. The construction of FEM solutions to elliptic interface 
problems dates back to 1970s [1^, and has been a subject of intensive investigation in 
the past few decades |T0l|3TJ [ITl • Since the elliptic interface problem deflned in Eqs. 
(1.1)- (1.6) provides opportunities to construct new FEM schemes, a wide variety of 
FEM approaches have been proposed in the literature. There are two major classes of 
FEM based interface methods, namely, interface-fltted FEMs and immersed FEMs, 
categorized according to the topological relation between discrete elements and the 
interface. In the interface-fltted FEMs, or body-fltted FEMs, the flnite element mesh 
is designed to align with the interface. Local mesh reflnement based priori and/or 
posteriori error estimation can be easily carried out. The performance of interface- 
fltted FEMs depends on the quality of the element mesh near the interface as well as 
the formulation of the problem. In fact, the construction of high quality FEM meshes 
for real world complex interface geometries is an active area of research. Bramble and 
King discussed a FEM for nonhomogeneous second order elliptic interface problems 
on smooth domains [4 . One way to deal with embedded interface conditions is to use 
distributed Lagrange multipliers [5], In fact, similar ideas have been widely used in 
mortar methods and flctitious domain methods for the treatment of embedded bound- 
aries [T3J. A Qi -nonconforming flnite element method was also proposed for elliptic 
interface problems [31j. Discontinuous Galerkin (DG) FEMs have been developed for 
elliptic equations with discontinuous coefficients [H [17] . Inherited from the original 
DG method, DG based interface schemes have the ffexibility to implement interface 
jump conditions. 

Immersed FEMs are also effective approaches for embedded interface problems 
[Tm \35\ \T5\ [21] ^ . A key feature of these approaches is that their element meshes 
are independent of the interface geometry, i.e., the interface usually cuts through 
elements. As such, there is no need to use the unstructured mesh to body-flt the in- 
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terface, and simple structured Cartesian meshes can be employed in immersed FEMs. 
Consequently, the time-consuming meshing process is bypassed in immersed FEMs. 
However, to deal with complex interface geometries, it is necessary to design appro- 
priate interface algorithms, which is similar to finite-difference based elliptic interface 
methods. In fact, immersed FEMs can be regarded as the Galerkin formulations of 
finite difference based interface schemes. It is not surprised that key ideas of many 
immersed FEMs actually come from the corresponding finite-difference based inter- 
face schemes. Additionally, other ideas in numerical analysis have been utilized for 
the construction of immersed FEMs. For example, a stabilized Lagrange multiplier 
method based on Nitsche technique has been used to enforce interface jump con- 
straints [16]. The performance of immersed FEMs depends on the design of elegant 
interface schemes for complex interface geometries. 

Convergence analysis of FEM based elliptic interface methods has been consid- 
ered by many researchers. Unlike the collocation formulation, the Galerkin framework 
of FEMs allows more rigorous and robust convergence analysis. Cai el al gave a proof 
of convergence for a DG FEM for interface problems [6 . Dryja et al discussed the 
convergence of the DG discretization of Dirichlet problems for second-order elliptic 
equations with discontinuous coefficients in 2D [9 . Hiptmair et al presented a conver- 
gence analysis of H{div] r^)-elliptic interface problems in general 3D Lipschitz domains 
with smooth material interfaces [19 . Recently, an edge-based anisotropic mesh re- 
finement algorithm has been analyzed and applied to elliptic interface problems [33 . 

Despite of numerous advancements in the numerical solution of interface prob- 
lems, there are still a few remaining challenges in the field. One of these challenges 
concerns the construction of higher order interface schemes. Currently, most interface 
schemes are designed to be of second order convergence. However, higher order meth- 
ods are efficient and desirable for problems associated with high frequency waves, such 
as electromagnetic and acoustic wave propagation and scattering, vibration analysis 
of engineering structures, and shock-vortex interaction in compressible ffuid ffows. It 
is easy to construct high order methods and even spectral methods for these prob- 
lems with straight interfaces in simple domains. However, it is extremely difficult to 
obtain high order convergence when the interface geometries are arbitrarily complex. 
A fourth order MIB scheme has been developed for the Helmholtz equation in media 
with arbitrarily curved interfaces in 2D [40 . Up to sixth order MIB schemes have been 
constructed for the Poisson equation with ellipsoidal interfaces in 3D |37j. There are 
two standing open problems concerning high order elliptic interface schemes, i.e., the 
construction of fourth-order 3D interface schemes for arbitrarily complex interfaces 
with sharp geometric singularities and the construction of sixth-order 3D interface 
schemes for arbitrarily curved smooth interfaces |37] . 

Another challenging issue in elliptic interface problems arises from nonsmooth 
interfaces or interfaces with Lipschitz continuity [39l [38l [36l [21] . Nonsmooth interfaces 
are also referred to as geometric singularities, such as sharp edges, cusps and tips, 
which commonly occur in real- world applications. It is a challenge to design high order 
interface schemes for geometric singularities both numerically and analytically. The 
first known second order accurate scheme for nonsmooth interfaces was constructed 
in 2007 [38]. Since then, many other interesting second order schemes have been 
constructed for this class of problems in 2D [8] [211 E] and 3D [37] • The second 
order MIB method for 3D elliptic PDEs with arbitrarily non-smooth interfaces or 
geometric singularities has found its success in protein electrostatic analysis [39 1 IT2 1 I7]. 
However, it appears truly challenging to develop fourth order schemes for arbitrarily 



4 



shaped nonsmooth interfaces in 3D domains, although fourth order schemes have 
been reported for a few special interface geometries [37] . Due to the need in practical 
applications, further effort in this direction is expected. 

The above mentioned difficulties in the solution of interface problems with geo- 
metric singularities significantly deteriorate in certain physical situations. It is well- 
known that the electric field diverges near the geometric singularities, such as tips of 
electrodes, antenna and elliptic cones, and sharp edges of planar conductors. Since 
electric field is related to the gradient of the electrostatic potential, i.e., the solution 
of the Poisson equation in the electrostatic analysis, it turns out that the solution of 
the elliptic equation has a lower regularity, e.g., the gradient does not exist at the 
geometric singularities. For isolated singularities, one can alleviate the difficulty by 
introducing an algebraic factor to solve a regularized equation whose solution has a 
high regularity [39, 36 . In this manner, second order MIB schemes have been con- 
structed in the past [39l[36]. However, it is also desirable to directly solve the original 
PDEs with a low solution regularity. The FEM developed by Hou et al is of first 
order convergence in the solution and 0.7th order convergence in the gradient of the 
solution when the solution of the Poisson equation is or continuous and the 
interface is Lipschitz continuous [21 . In fact, due to the low solution regularity in- 
duced by the nonsmooth interface, it is very difficult to analyze the convergence of 
numerical schemes because commonly used techniques may no longer be available. 
Therefore, there is a pressing need to develop both high order numerical methods 
and rigorous analysis of numerical methods for elliptic interface problems with low 
solution regularities induced by geometric singularities. 

The objective of the present work is to construct a new numerical method for ellip- 
tic equations with low solution regularities induced by nonsmooth interfaces. To this 
end, we employ a newly developed weak Galerkin finite element method (WG-FEM) 
by Wang and Ye [34 . Like the discontinuous Galerkin (DG) methods, WG-FEM 
makes use of discontinuous functions in the finite element procedure which endows 
WG-FEMs with high flexibility to deal with geometric complexities and boundary 
conditions. For interface problems, such a flexibility gives rise to robustness in the 
enforcement of interface jump conditions. Unlike DG methods, WG-FEM enforces 
only weak continuity of variables naturally through well defined discrete differential 
operators. Therefore, weak Galerkin methods avoid pending parameters resulted from 
the excessive flexibility given to individual elements. As a consequence, WG-FEMs are 
absolutely stable once properly constructed. Recently, WG-FEMs have been applied 
to the solution of second-order elliptic equations [26^ , and the solution of Helmholtz 
equations with large wave numbers [27]. A major advantage of the present weak 
formulation is that it naturally enables WG-FEM to handle interface problems with 
low solution regularities. We demonstrate that the present WG-FEM of lowest order 
is able to achieve from 1.75th order to second order of convergence in the solution 
and about first order of convergence in the gradient when the solution of the elliptic 
equation is or H'^ continuous and the interface is or Lipschitz continuous. 

The rest of this paper is organized as follows. Section |2] is devoted to a description 
of the method and algorithm. We shall design a weak Galerkin formulation for the 
elliptic interface problem given in Eqs. ([lT])-([L6|. Section [s] is devoted to a con- 
vergence analysis for the weak Galerkin scheme presented in Section [2j In Section [4] 
we shall present some numerical results for several test cases in order to demonstrate 
the performance of the proposed WG-FEM for elliptic interface problems in 2D. We 
first consider a few popular benchmark test examples with smooth but complex in- 
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terfaces. We then carry out some investigation about nonsmooth interfaces. Some 
of these test problems admit solutions with low regularities, for which our numerical 
results significantly improve the best known result in literature and are better than 
our theoretical prediction. This paper ends with a conclusion, especially a remark on 
the use and possible advantages of high order WG-FEMs. 

2. Methods and Algorithms. Let 7/i be a partition of the domain Q with 
mesh size h. We require that the edges of the elements in Th align with the interface 
r. Thus, the partition Th can be grouped into two sets of elements denoted by 
Tfj; = 7/1 n and = 7/i fl (^2, respectively. Observe that 7^ provides a finite 
element partition for the subdomain flj^j = 1, 2. The intersection of the partition Th 
also introduce a finite element partition for the interface F, which shall be denoted 
by Th' For simplicity, we adopt the following notation: 



{v^w)k '= / vwdK^ {v^w)dK = / vwds. 

Jk JdK 

For each triangle G 7/i, let and dK denote the interior and boundary of 
K respectively. Denote by Pj{K^) the set of polynomials in with degree no more 
than j, and P£{e) the set of polynomials on each segment (edge or face) e, e G dK 
with degree no more than £. A discrete function w = {wq, wi^} refers to a polynomial 
with two components in which the first component wq is associated with the interior 
Kq and wi) is defined on each edge or face e, e G dK. Please note that Wb may or may 
not equal to wq on dK. Now we introduce three trial finite element spaces as follows 

(2.1) Uh -{w = {wo,Wb} : {wo,Wb}\K e Pj{K°) x Pe{e), e € dK,\/K e T,^} , 

(2.2) Vh :={p = {po,Pb} : {po,Pb}\K € Pj{K°) x P^(e), e G dK,yK G T^} , 

(2.3) A,, --{p : pU e -Pm(e), e G . 

Define two test spaces by 

Uh = {w = {wo, Wb} e Uh ■■ Wb\e = 0, for e G \ T} 

and 

Vh ={P= {Po, Pb} eVh-. Pb\e = 0, for eed^2\ T}. 

For each w = {k;o,'^6} ^ Uh or Vh^ we define the discrete gradient of denoted by 
V dvo G Vr{K) on each element K^ by the following equation: 

(2.4) / VdW • qdK = - wo{V • q)dK + / wt{q • n)ds, \fq G Vr{K), 
Jk Jk JdK 

where Vr{K) is a subspace of the set of vector- valued polynomials of degree no more 
than r on i^. 

The selection of the indices j, m, and r is critical in the design of weak Galerkin 
finite element methods. A detailed discussion on the selection of those indices can be 
found in [34 . In the present study of the interface problem, we shall let j = ^ = m = 



/c > in (2.1)-(2.3) and chose the Raviart-Thomas element for Vr{K) := RTj.{K). 
These elements are referred as {Pfe(i^°)^ 7^fe(e)^ ^fe(r)} element in our numerical 
experiments. An exploration of other possible combinations is left to interested readers 
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for future research. Recall that the Raviart- Thomas element [32] RTk{K) of order k 
is of the following form 

RTk{K) = Pk{Kf + Pk{K)x, 

where Pk{K) is the set of homogeneous polynomials of degree k and x = (xi,X2). 

Weak Galerkin Algorithm 1. A numerical approximation for 
can be obtained by seeking Uh = {uq^ui)} G Uh satisfying = QbSi on dQi \ T, 
= {^0^ Vb} eVh satisfying = Qbg2 on dVt2 \ T and G such that 

(2.5) {AVdUh. Vdw) - {Xh, Wb)r= (/i, ^o), G 

(2.6) {AWdVh. VdP) + (A/,, P6)r= (/2, Po) + (V^, />5)r, Vp G F;? 

(2.7) {ub - Vb, /i)r= ((/>, M)r, V/i G A^. 

i/ere Qbgi ctnd Qbg2 ct'i^e the standard projection of the Dirichlet boundary data in 
Pk{e) for any edge/ face e G dVt. 

3. Convergence Theory. The goal of this section is to provide a convergence 
theory for the weak Galerkin finite element method as described in the previous 
section. First, we show that the WG algorithm has one and only one solution in the 
corresponding finite element trial space. 

Lemma 3.1. The weak Galerkin finite element method (2.5)-(2.1) has a unique 
solution. 

Proof. It suffices to show that the solution of (2.5)-(2.7) is trivial when data is 
homogeneous; i.e., /i = /2 = = ^2 = ^ = = 0. To this end, by letting /j, = Ub — Vb 
in (2.7) and then using the assumption of = we obtain Ub = Vb on T. Next, by 
first letting w = Uh in (2.5) and p = Vh in (2.6) and then adding them up gives 

{AVdUh, VdUh) + {AVdVh^ VdVh) = 0, 

which implies S/dUh = and S/dVh = 0- Thus, both and are constants over Qi 
and O2, respectively. It follows from the assumption of Ub = on the non-empty set 
dQi \ r that Uh = {uo^Ub} = 0. Thus, we have Vb = Ub = on F, which together 
with \/dVh = implies Vh = 0. Finally, using VdUh = and letting Wb = Xh on any 
edge/face e G F/^,, we have from equation (|2.5|) that 



(A/i, Xh)r = 0, 

which implies Xh = 0. □ 

Denote by Qh = {Qo^Qb} a local projection operator where: 

Qo : H\K^) ^ Pk{K^), Qb : H^e) ^ Pfc(e), e G dK 

are the usual projections into the corresponding spaces. Let 11^ be the usual 
projection operator in the mixed finite element method such that Il^q G H^div^Qi); 
and on each K eTh^ one has Il/^q G RTk{K) satisfying 

(V • q, Wo)k = (V • n/,q, Wo)k. V^o e Pfe(i^°). 

Lemma 3.2. Let r G i^(div, 1]^) be a smooth vector-valued function and Uh be the 
locally defined projection operator commonly used in the mixed finite element method. 
Then, the following identify holds true 

(3.1) (-V • r, wo)k = ^ (n/,r, Vd^)^ - (^6, r • n)r 

Ker' Ken 
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for any w = {wq, w^} G V with V = or and z = 1, 2. 
Proof. It follows from the definition of Hh and \/d that 

^ (-V • r, wo)k = ^ (-V • UhT, wo)k 

Ken 

Now using the continuity of HhT ■ n across each interior edge or face we arrive at 
^ (-V • r, wo)k = Vd^)^ - {wb, UhT • n)QQ^ 

Ken 

= (n^r, \/dw)K - {wb, T ■ n)r, 
Ken 

where we have used the fact that {w^^ H^r • n)r = {w^^ t • n)r and i^;^ = on dfti\T. 
This completes the proof. □ 

The following error estimates are straightforward from the definition of Hh and 
Qh- Readers can also find a verification of the result in [34J. 

Lemma 3.3. For u e H^^'^(Vti) and v e H^^'^(yt2) with k>0, we have 

(3.2) \\Ilh{A\/u) - A\/d{Qhu)\\ < C/i^+iix|U+2,Q,, 

(3.3) \\Uh{A\/v) - A\/d{Qhv)\\ < Ch'^'\\v\\k+2,n,^ 



It is well known that there exists a constant C such that for any function g G 



(3.4) 



\9\\l<C{h-'\\g\\l^hK\\g\\i 



K) : 



where e is an edge of K. 

Theorem 3.4. Let {uh^Vh^Xh) ^ Uh x Vh x Ah be the solution arising from the 
weak Galerkin finite element scheme (2.5)-(2.1). Then, the following error estimates 
hold true 

(3.5) \\Wd{QhU - Uh)\\ + \\Vd{QhV - Vh)\\< Ch''+\\\u\\k+2 + \\v\\k+2), 

(3.6) \\AWu ■ n - A„||r< Ch''+'^\\u\\k+2 + MM, 

where C stands for a generic constant independent of the mesh size h. 

Proof For simplicity, we assume that the coefficient A is a constant tensor on 
each element K. The proof can be extended to the general case of non-constant tensor 
A without any difficulty. 

Now testing (1.1) with w e and then using (3.1) leads to 

(UhiAVu), Vdw) - {AVu ■ n, ^5)r = (/i, ^o). 
Adding the term {A\/dQhU^ V^k;) to both sides of the equation above, we obtain 

{AVdQhU, Vd^) - {A\/u • n, ^6)r= (/i, ^o) 

(3.7) -{Uh{AVu) - AVdQhU, Vdw). 



Similarly, on Q2 we have for any p G the following relation 



{AVdQhV, Vdp) - {AVv ■ n, pb)r= (/2, Po) 

(3.8) -{Uh{A\/v) - AWdQhV, VdP)- 

Using the interface condition ( |1.6| ), we from from ( |3.8| ) that 

{AVdQhV, VdP) + {AVu • n, p5)r= (/2, Po) + (V^, P6)r 

(3.9) -(n;,(AV^) - AVdQhV. ^dP)^ 



Testing the interface condition (1.5) by p G Kh implies 

(3.10) {QhU - QbV, p)r = {u-v, p)r = /i)r. 

For simplicity, we introduce the following notations to represent the errors: 



= {^05 ^b} 



= {Qou - uo, QhU - iib}, 
= Qb{A\/u ■ n) - Xh. 



Thus, the difference of ( |3.7[ ) and ( |2.5[ ) gives 

(3.11) (AVrfe)^, Vrf^) - (e,, ^5)r = -{Uh{AVu) - AVdQhU, Vdw). 
The difference of (|3^ and ([2^ gives 

(3.12) {AVdel VdP) + (en, P6)r = -(n;,(AV^) - AVdQhV. "^dP)^ 
The difference of ( |3.1Q[ ) and ([2^71 gives 

(3.13) (e^-e^, M)r = 0, V/i G A/,. 



First letting w = e'^ and p = e]^ in (3.11 ) and (3.12) and then adding them up gives 

||A^V,e;^|p + P^V,ej;|p = -{Uh{AVu)-AVdQhU, V del) 
-{Iih{AVv)-AVdQhV, Vdel). 



Using ( |3.2[ )-( |3.3| ), we have from the above equation that 

C{\\Vael\\ + \\yael\\f< ||A^V,e^f + V.e^f 

< \\Ilh{AVu)-AVdQhu\\ WVdelWW 
+ ||n^(^V^;)-^VdQ^^;|| WV^elW 

< Ch'^+\\\uU+2 + ||f IU+2)(||V,e^|| + ||V,e;;||), 



which implies 
(3.14) 



W^delW + WVdelW < Ch'^+H\\u\\k+2 + \\vh+2). 



Using (3.4) and the inverse inequality, we have for any edge or face e C F fi dK 



(3.15) 
(3.16) 



\\AVdel-n\\l<Ch-'\\Vdel\\l, 
(n/,(^Vw) - AVdQhu) ■ n\\l< Ch-^\\n,^{AVu) - AVdQhu\\%. 
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Choosing w = {wo^Wb} in equation (3.11) such that wq = for ah K e Tfj; and 
W}) = for e G r and w^j = otherwise yields 



Using (2.4), the equation above becomes 



Iknllr = ^{{A^dcl • n, eu)e + {{Uh{AVu) - AVdQhu) ■ n, e^)e). 



eGr 



Using ( [3l5l ), ( [3l6l ), ( [3^ and ([3^, we have 



|BC^(||AV,e^n|| 



eer 



{Uh{A\/u) - AWdQhu) • n||e||e,||e) 

1/2 



hU\\K 



en r 



l^ll/c+2 + ||^||/c+2j||en||r, 



which verifies the error estimate (3.6). This completes the proof of the theorem. □ 



4. Numerical experiments. The goal of this section is to numerically validate 
the proposed WG algorithm by solving some benchmark elliptic interface problems 
for which analytical solutions are known. To fully demonstrate the accuracy and 
robustness of the WG method, we consider challenging problems involving Lipschitz 
continuous interfaces, highly oscillated solutions, and solutions with low regularities. 
For simplicity, we use the piecewise constant finite element Po{dK) — Pq{K) — RTo{K) 
on structured triangular meshes in all test cases. The mesh generation and computa- 
tion are all conducted in the MATLAB environment. Since the analytical solutions 
are known for each test case, the nonhomogeneous terms /i and /2 of the elliptic equa- 
tions and the Dirichlet boundary data can be correspondingly derived. Moreover, the 
two interface jump conditions for the solution and the fiux across the interface are 
calculated according to the given analytical solution. Numerical errors in the solution 
and its gradient are reported in Loo norms in all examples. 



Table 4.1 

Numerical convergence test for Example 1. 



Mesh 


max{/i} 


Solution Gradient 


Loo error order 


Loo error order 


Level 1 


2.8553e-01 


1.0266e-02 


1.3042e-02 


Level 2 


1.5110e-01 


2.9631e-03 1.9525 


6.5234e-03 1.0886 


Level 3 


7.7543e-02 


7.6763e-04 2.0247 


3.2615e-03 1.0391 


Level 4 


3.9258e-02 


1.9518e-04 2.0118 


1.6306e-03 1.0185 


Level 5 


1.9749e-02 


4.9198e-05 2.0058 


8.1523e-04 1.0090 



Example 1. We first study a classical circular interface problem [42 . Consider 
a square domain [—1,1] x [—1,1] with a circular interface = ^ = |. The 
coefficient A is defined to he Ai = b and ^2 = 2, respectively on each subdomain, for 
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-1 -1 




-1 -1 



Fig. 4.1. WG solutions of Example 1. Left: Mesh level 1; Right: mesh level 5. 



r > 0.5 and r < 0.5. The analytical solution to the elliptic equation is given as 
"1 



1 1 

4 V' ~ Sb~ b 
v{x,y) = -{x^ - 1) 



r > 0.5 
r < 0.5. 



By choosing b = 10, the function and flux jumps are actually constants [42\ 

By using a uniform triangular mesh, the Loo error for the solution and its gradient 
of the WG method is reported in Table 4.1 Based on successive mesh reflnements. 



the numerically detected convergence rates are also reported for both error measure- 
ments. It can be seen that the orders of convergence in Loo norm for the solution 
and gradient are, respectively, two and one for piecewise constant WG flnite element 
approximations. This agrees with our theoretical results. 

In Fig. |4.1[ the WG solutions based on mesh levels 1 and 5 are depicted. By using 
piecewise constant flnite elements, the WG solution at mesh level 1 clearly consists of 
piecewise constants. On the other hand, a much better numerical solution is attained 
at mesh level 5. It is also seen that the function and flux jumps are constant across 
the circular interface. 

Table 4.2 

Numerical convergence test for Example 2 with k, = 2. 



Mesh 


max{/i} 


Solution Gradient 


Loo error order 


Loo error order 


Level 1 


3.1778e-01 


7.6628e-02 


1.5852e-01 


Level 2 


1.5889e-01 


1.5119e-02 2.3415 


5.2258e-02 1.6009 


Level 3 


7.9444e-02 


3.6142e-03 2.0646 


2.0731e-02 1.3338 


Level 4 


3.9722e-02 


8.9375e-04 2.0157 


9.6505e-03 1.1031 


Level 5 


1.9861e-02 


2.2492e-04 1.9905 


4.7340e-03 1.0275 



Example 2. To further explore the potential of the proposed WG method, 
we consider a circular interface problem with highly oscillatory solutions. The high 
frequency or short wave solutions are commonly encountered in solving electromag- 
netic wave propagation and scattering problems governed by time domain Maxwell's 
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equations or frequency domain Helmholtz equations. Here, we consider the following 
Helmholtz interface problem |42] 



(4.1) 
(4.2) 



-V • (AiS/u) ■ 



klv 



h 

/2 



r > 0.5, 
r < 0.5, 



where the interface is still r'^ = \ over the domain [—1, 1] x [—1, 1]. The wavenumber 
ki = ncFi {i = 1,2) depends on the wavenumber of the wave solution in free space 
and dielectric profiles cr^. In the present example, both A and are chosen to be 
discontinuous across the interface. In particular Ai = 1 and ai = \/T0 for r > 0.5 
and ^2 = 10 and (J2 = 1 for r < 0.5. The analytical solution is given as 



uix^y) = —sm{f<ix) cos(/^y) 



r > 0.5, 
r < 0.5. 



Other necessary interface jump conditions can be derived from the analytical solution. 

The WG solution of the Helmholtz equation with high wavenumbers has been 
explored in [27 . For the present Helmholtz interface problem, a WG algorithm can be 
similarly constructed based on equations (2.5) - (2.7). Two frequencies are tested, i.e., 
K = 2 and k 



the solution is of high frequency, see Fig. 4.2 



8, and the results are reported in Tables [42] and [431 using k, = 

In comparison with the numerical 
accuracy of the WG method for both cases, the WG errors for = 8 are larger 
even though a finer mesh is employed in the same mesh level. Nevertheless, when we 



Table 4.3 

Numerical convergence test for Example 2 with = 8. 



Mesh 


max{/i} 


Solution Gradient 


Loo error order 


Loo error order 


Level 1 


1.5837e-01 


8.6774e-01 


5.2751e+00 


Level 2 


7.9186e-02 


1.4156e-01 2.6159 


9.1695e-01 2.5243 


Level 3 


3.9593e-02 


2.3768e-02 2.5743 


2.8679e-01 1.6768 


Level 4 


1.9796e-02 


5.8719e-03 2.0170 


1.0732e-01 1.4180 


Level 5 


9.8982e-03 


1.4717e-03 1.9964 


4.7255e-02 1.1834 





I 



-1 -1 



Fig. 4.2. WG solutions at mesh level 5 for Example 2. Left: k = 2; Right: = 8. 
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examine the convergence rates, the WG method maintains the same rates, i.e., second 
and first orders, respectively, for the solution and its gradient, in both cases. This 
verifies the capability of the WG method in resolving short wave solutions. 



Table 4.4 

Numerical convergence test for Example 3 with 6 = 10. 



Mesh 


max{/i} 


Solution Gradient 


Loo error order 


Loo error order 


Level 1 


5.6522e-01 


1.4339e-01 


1.5956e-01 


Level 2 


2.8261e-01 


2.6979e-02 2.4100 


9.1557e-02 0.8014 


Level 3 


1.4130e-01 


7.3332e-03 1.8792 


5.8167e-02 0.6544 


Level 4 


7.0652e-02 


1.4904e-03 2.2988 


2.9623e-02 0.9735 


Level 5 


3.5326e-02 


2.8124e-04 2.4058 


1.5250e-02 0.9579 




-1 -1 -1 -1 



Fig. 4.3. WG solutions of Example 3 with h = 10. Left: Mesh level 1; Right: mesh level 5. 



Table 4.5 

Numerical convergence test for Example 3 with h = 1000. 



Mesh 


max{/i} 


Solution Gradient 


Loo error order 


Loo error order 


Level 1 


5.6522e-01 


4.1482e-01 


3.1081e-01 


Level 2 


2.8261e-01 


5.3753e-02 2.9481 


1.0784e-01 1.5271 


Level 3 


1.4130e-01 


1.3707e-02 1.9713 


5.6088e-02 0.9431 


Level 4 


7.0652e-02 


3.0677e-03 2.1598 


2.9658e-02 0.9193 


Level 5 


3.5326e-02 


6.2057e-04 2.3055 


1.5204e-02 0.9640 



Example 3. We next consider an example with a high contrast in coefficient A 
[42] to examine the robustness of the WG method. Consider the second order elliptic 
equation over the domain [—1, 1] x [—1, 1]. The interface F is defined to be an ellipse 
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We set A2 = b inside of F, and Ai = 1 outside of T. The analytical solution is given 
as 



u{x,y) = 5e ^ ^ 
v{x^ y) = cos(y) 



outside r, 
inside F. 



Other necessary conditions can be derived from the analytical solution. 

Two values of b are tested with b = 10 and b = 1000. The latter case involves a 
much larger jump in A. This results in a larger condition number for the corresponding 
discretized linear system. This means that extensive computational time is usually 
needed in the matrix solving. The WG method performs robustly for both cases. The 
numerical results are presented in Table 4.4 and Table |4.5[ respectively, for b = 10 
and b = 1000. It can be seen in both tables that the convergence rates are first and 
second orders, respectively, for the gradient and the solution in Loo norm. The WG 



solutions on mesh levels 1 and 5 are shown in Fig. |4.3| and Fig. |4.4[ respectively, for 
b = 10 and b = 1000. We note that in both cases, the analytical solutions are the 
same and are independent of the coefficient A. Thus, even though there are some 
minor differences in the WG solutions at mesh level 1 for 6 = 10 and b = 1000, the 
WG solutions are visually identical at mesh level 5. This demonstrates the robustness 
of the WG method in handling interface problems with high contrast. 

Table 4.6 

Numerical convergence test for Example 4- 



Mesh 


max{/i} 


Solution Gradient 


Loo error order 


Loo error order 


Level 1 


3.4726e-01 


1.5885e-03 


3.5260e-02 


Level 2 


1.7363e-01 


3.8053e-04 2.0616 


1.6630e-02 1.0842 


Level 3 


8.9010e-02 


9.1082e-05 2.1399 


9.0379e-03 0.9126 


Level 4 


4.9602e-02 


2.6778e-05 2.0936 


5.0265e-03 1.0034 


Level 5 


2.4286e-02 


5.9847e-06 2.0982 


2.4459e-03 1.0087 



Example 4. We next consider a classical elliptic interface problem with both 
concave and convex curve segments [42 . The interface F is parametrized with the 
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Fig. 4.4. WG solutions of Example 3 with h = 1000. Left: Mesh level 1; Right: mesh level 5. 
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polar angle 

_ 1 sin(5l9) 

inside a square domain [—1, 1] x [—1, 1]. The coefficient is chosen to be Ai = 10 and 
^2 = 1, respectively, for outside and inside F. The analytical solution is given as 

u{x, y) = 0.1(x^ + y^)'^ - 0.01 ln(2 V^^T^) outside T 

2, 2 

v{x^y) = inside F. 

Other necessary interface jump conditions can be derived from the analytical solutions. 



The WG solutions and their convergence are reported in Table [46] and Fig. |4.5 
respectively. It can be seen from the table that the convergence is of order two and one 
for the solution and its gradient, respectively, in the L^q norm. It should be pointed 
out that mesh quality plays an important role in the performance of general finite 
element methods. For weak Galerkin, we noted that the quality of mesh affects the 
corresponding numerical approximation in Loo norms. For example, we noted a slower 
convergence when the finite element partition consists of both acute and non-acute 
triangles with nontrivial portion. After eliminating most of the non-acute triangles, 
the WG approximation exhibits an improved rate of convergence numerically. Figure 



4.6 demonstrates some meshes generated by MATLAB in the present study after an 
elimination of most non-acute triangles through a numerical procedure. It can be seen 
that, at each level, the triangular mesh is fitted to the interface at the nodal points. 
Observe that these meshes are highly non-uniform and unstructured. In particular, it 
is clear that much smaller triangles are employed to resolve the concave portion of the 
interface. This somehow helped the WG method to deliver numerical solutions with 
high accuracy since the solution changes rapidly in regions with large curvatures. 



Example 5. In this and next two examples, we consider an interface problem 
with continuous interfaces as discussed in [21j. The domain is given as 1] = 
[—1,1] X [—1, 1], and the interface F is defined with two pieces: y = 2x ioi x ^ y > {) 
and 7/ = 2x + x^, for x + y < 0. The interface F divides the domain into two parts 
Vti and 0.2' Here we denote Vti and 1^2, respectively, to be the part on the left and 




Fig. 4.5. WG solutions of Example 4- Left: Mesh level 1; Right: mesh level 5. 
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right of r. An inhomogeneous second order elliptic equation is considered, i.e., the 
coefficient A is not piecewise constant. Instead, we have Ai{x^y) = {xy + 2)/5 in Qi 
and A2{x^y) = (x^ — y'^ -\- 3)/7 in r^2- The analytical solution u is fixed to be = 2 




Fig. 4.6. Finite element meshes of Example 4- Top left: Mesh level 1; Top right: Mesh level 
3; Bottom: Mesh level 5. 
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in 1^1, while in 1^2, ^ is defined piecewisely by 
v{x,y) 



sm{x + y), if X -\- y < 
X -\-y^ if X -\- y > 0. 



Other necessary conditions can be derived from the analytical form of the solution. 
In particular, the Dirichlet boundary data is derived for both u = g and v = since 
in this example, Q2 is not enclosed by Qi. 



A numerical investigation was conducted and the results are reported in Table [477 
Again, the WG method attains the theoretical order of convergence for this example, 
i.e., first and second order of accuracy respectively for the gradient and the solution 
in Loo norm. We note that due to the piecewise definition, the solution v{x^y) is 

continuous but not across the line x -\- y = 1. In our computation, the line 
X -\- y = 1 is not treated as a boundary. Thus, our triangulation is not aligned with 
this line. This can be seen clearly in the graph of the WG solution on mesh level 1 



(Fig. 4.7 left). In the right chart of Fig. 4.7, a smooth transition can be seen across 



X -\- y = 1 in the WG solution on mesh level 5. 



Example 6. Here we consider the same interface and domain geometry as in 
Example 5. The coefficient function A{x^y) is defined in the same way as in Example 
5. We also fix the solution as = 2 in l^i. But in (^2, the solution is given differently 



Table 4.7 

Numerical convergence test for Example 5. 



Mesh 


max{/i} 


Solution Gradient 


Loo error order 


Loo error order 


Level 1 


6.2575e-01 


L5275e-02 


4.8833e-02 


Level 2 


2.9662e-01 


3.9259e-03 1.8200 


3.1059e-02 0.6062 


Level 3 


L5170e-01 


L0586e-03 L9546 


1.5584e-02 1.0285 


Level 4 


7.9102e-02 


3.3419e-04 L7707 


8.0209e-03 1.0200 


Level 5 


4.4180e-02 


1.1698e-04 1.8022 


4.4282e-03 1.0199 




Fig. 4.7. WG solutions of Example 5. Left: Mesh level 1; Right: mesh level 5. 
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by (see [21] for details) 



v{x,y) 



sm{x -\- y) -\- cos{x + ?/), 
x + + 1, 



if X -\- y < 
if X + ^ > 0. 



Other necessary conditions can be derived similarly. Observe that the solution v = 
v{x^ y) is continuous but not across the line x + ^ = 1. 

The WG algorithm is formulated as in the previous example. In particular, since 
v{x^y) is across x -\- y = 1, the function and flux jumps are trivially zero across 
x-\-y = 1. Thus, again, no special boundary or interface treatment is carried out near 
the line x -\-y = 1. In fact, this line actually cuts through the finite element triangles; 



see the left chart of Fig. 4.8 It can be seen that the loss of regularity in v{x,y) is 
visually indistinguishable in the WG solution on mesh level 5 (see the right chart of 
Fig. 4.8). In literature, a certain reduction of convergence order has been reported for 
this continuous example [21 when the standard Galerkin finite element method 
is employed for a numerical approximation. In the present study, it is found that the 
order of convergence remains unchanged for the WG method, see Table 4.8 This 
indicates that the WG method preforms robustly for the challenging problem with 
low regularity. 

Example 7. Consider again the same interface problem as in Examples 5 and 
6, over a larger domain Q = [—1,3] x [—1, 1]. The coefficient function A{x^y) is now 
defined to be Ai{x^ y) = 1 in and A2{x^ y) = 2 -\- sin{x + y) in • The analytical 



Table 4.8 
Numerical convergence test for Example 



Mesh 


max{/i} 


Solution Gradient 


Loo error order 


Loo error order 


Level 1 


6.2575e-01 


2.6559e-02 


4.9273e-02 


Level 2 


2.9662e-01 


5.9674e-03 2.0001 


2.8456e-02 0.7355 


Level 3 


1.5170e-01 


1.5342e-03 2.0257 


1.4493e-02 1.0062 


Level 4 


7.9102e-02 


4.4127e-04 1.9137 


7.5081e-03 1.0100 


Level 5 


4.4180e-02 


1.5022e-04 1.8500 


3.9584e-03 1.0990 
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-1 -1 
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Fig. 4.8. WG solutions of Example 6. Left: Mesh level 1; Right: mesh level 5. 
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solution is give as [21] 

u{x^ y) = S in r^i, 

v{x, y) = {x^ + y'^f^^ + sin(x + y) in 1^2. 

Other necessary conditions can be derived from the analytical solutions. 

The analytical solution is piecewise in this example. In particular, the function 
v{x^ y) has a singularity at (0, 0) with blow-up derivatives. Even though the analytical 
solution is well behaved in other places, the singularity at the origin may cause some 
order reduction in the finite element simulation [21 . Due to this singularity, the 
convergence rate of the WG method in the Loo norm of the solution is reduced to 
be about 1.75 for the present example, see Table [43) On the other hand, the order 
of convergence for the WG method for the gradient remains to be about one. The 



numerical solutions at mesh levels 1 and 5 are depicted in Fig. 4.9 Interestingly, in 
both levels, the singularity at the origin is invisible. But the numerical algorithm can 
sense this singularity via the convergence rate. 

Example 8. In this and next two examples, we continue to examine the perfor- 
mance of the WG method for solutions of low regularity when the interface is fixed to 
be a Lipschitz continuous curve [ST. On a rectangular domain Vt = [—1, 3] x [—1, 1], we 
consider a case where the interface T consists of two pieces: y = 2x ioi x ^ y > {) and 
y = — x/2, for X -\- y < 0. Thus, the interface T has a kink at (0, 0). Denote by Qi the 
region on the left and upper part of F, and Q2 to be the rest of the domain. The coeffi- 



Table 4.9 

Numerical convergence test for Example 7. 



Mesh 


max{/i} 


Solution Gradient 


Loo error order 


Loo error order 


Fevel 1 


6.5801e-01 


3.4340e-02 


3.2817e-01 


Fevel 2 


3.4551e-01 


1.1022e-02 1.7641 


1.7896e-01 0.9413 


Fevel 3 


1.7932e-01 


3.5197e-03 1.7405 


9.5039e-02 0.9650 


Fevel 4 


9.1176e-02 


1.0806e-03 1.7459 


5.0095e-02 0.9468 


Fevel 5 


4.6365e-02 


3.3092e-04 1.7499 


2.6003e-02 0.9696 




Fig. 4.9. WG solutions of Example 7. Left: Mesh level 1; Right: mesh level 5. 
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cient is given piecewisely hy Ai{x^y) = {xy-\-2)/5 in l^i and A2{x^y) = (x^— ?/^ + 3)/7 
in The analytical solution u is fixed to be = 8 in (^i, while in (^2, is given 
piecewisely by 



v{x,y) 



s'm{x -\- y)j if X -\- y < 0^ 
X -\- y^ if X -\- y > 0. 



Other necessary conditions can be derived from the analytical solutions. 

Similar to Example 5, the analytical solution is of continuous but not in 
this example. The WG method yields uniformly second and first order of accuracy in 
Loo norm, respectively, for the solution and its gradient as shown in Table 4.10 The 



kink of the interface F is clearly seen in both charts of Fig. 4.10 In particular, on 
the coarsest mesh, the alignment of finite element triangles to two line segments of F 
is obvious. No deterioration in solution occurs around the interface corner. 

Example 9. Consider an elliptic problem with the same interface and domain 
geometry as in Example 8. The coefficient function A{x,y) is defined the same as in 
Example 8. Also, we fix the solution to be = 8 in l^i. But in 1^2, the solution is 
given differently [21^ by 



v{x,y) 



sin{x -\- y) -\- cos{x + ?/), 
x + + 1, 



if X -\- y < 0, 
if X + 7/ > 0. 



Table 4.10 
Numerical convergence test for Example 8. 



Mesh 


max{/i} 


Solution Gradient 


Loo error order 


Loo error order 


Level 1 


5.6919e-01 


1.4984e-02 


3.3001e-02 


Level 2 


2.8459e-01 


4.4307e-03 1.7578 


1.7812e-02 0.8897 


Level 3 


1.4230e-01 


1.2336e-03 1.7687 


9.8402e-03 0.8561 


Level 4 


7.1149e-02 


3.2639e-04 1.8634 


5.1905e-03 0.9228 


Level 5 


3.5574e-02 


8.5110e-05 1.9231 


2.6699e-03 0.9591 





I 



Fig. 4.10. WG solutions of Example 8. Left: Mesh level 1; Right: mesh level 5. 



20 



Other necessary conditions can be derived similarly. Like in Example 6, v{x,y) is of 
continuous but not across the line x -\- y = 1. Similarly, no special numerical 
treatment is invoked near the line x -\- y = 1 so that this line actually cuts through 



finite element triangles, see the left chart of Fig. 4.11 It can be seen from Table 4.11 
that the WG method converges uniformly with orders being close to second and first, 
respectively, for the solution and its gradient in the Loo norm. Technically speaking, 
the interface is Lipschitz continuous in this example, while it is continuous in 
Example 6. In other words, the interface regularity is lower here than that in Example 
6. However, since body-fitted triangular meshes are employed in the WG method, the 
interface with two simple line segments in the present example is easier to be resolved 
by the body- fitted grids than the curved interface in Example 6. The only issue that 
could make an impact on the convergence is the geometrical singularity at the origin. 
The last and present examples show that the WG method is very robust in handling 
geometrical singularity, or sharp edged corners in general. 

Example 10. Consider the same interface and geometry as in Examples 8 and 
9. The coefficient function A{x,y) is now defined to be Ai{x,y) = 1 in l^i and 
A2{x,y) = 2 -\- sm{x + y) in The analytical solutions are give as [21] 



u{x,y) 8 

v{x^ y) = (x^ + y^)^^^ + sin(x + y) 



in 1^1, 
in 0.2- 



Other necessary conditions can be derived from the analytical solution. 



Table 4.11 
Numerical convergence test for Example 9. 



Mesh 


max{/i} 


Solution Gradient 


Loo error order 


Loo error order 


Level 1 


5.6919e-01 


1.4326e-02 


3.3514e-02 


Level 2 


2.8459e-01 


4.4250e-03 1.6949 


2.0051e-02 0.7411 


Level 3 


1.4230e-01 


1.2350e-03 1.8411 


1.0768e-02 0.8968 


Level 4 


7.1149e-02 


3.2693e-04 1.9175 


5.5458e-03 0.9573 


Level 5 


3.5574e-02 


8.5118e-05 1.9415 


2.8004e-03 0.9858 




i 



Fig. 4.11. WG solutions of Example 9. Left: Mesh level 1; Right: mesh level 5. 
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Similar to Example 7, the analytical solution is piecewise in the present case. 
The solution v{x^y) has a singularity at (0,0) with blow-up derivatives. The WG 
method achieves a similar order of accuracy as that in Example 7. In particular, it 



can be seen in Table 4.12 that the L^o error in the solution is of order 1.75, while it 
remains to be about 1 for the gradient approximation. Both function singularity and 
geometrical singularity are presented at the origin in this test; see the WG solutions 
in Fig. |4.12[ Nevertheless, the WG method is capable of delivering decent results at 
mesh level 1. 

Finally, we would like to summarize the WG results from Example 5 to Example 
10 in Table |4.13| In each case, an overall convergence rate is reported, which is 
calculated based only on mesh level 1 and level 5. It can be seen from the Loo error 
of the gradient that the WG method essentially attains the first order of accuracy 
in all test cases. For the Loo error of the solution, we analyze the order of the WG 



Table 4.12 
Numerical convergence test for Example 10. 



Mesh 


max{/i} 


Solution Gradient 


Loo error order 


Loo error order 


Level 1 


7.1020e-01 


4.0630e-02 


7.1308e-02 


Level 2 


3.5510e-01 


1.1389e-02 1.8349 


4.0849e-02 0.8038 


Level 3 


1.7755e-01 


3.3569e-03 1.7624 


2.1170e-02 0.9483 


Level 4 


8.8775e-02 


1.0150e-03 1.7256 


1.0780e-02 0.9737 


Level 5 


4.4387e-02 


3.1128e-04 1.7052 


5.5740e-03 0.9516 




Fig. 4.12. WG solutions of Example 10. Left: Mesh level 1; Right: mesh level 5. 

Table 4.13 

Summary of the overall orders of the WG method for solutions with low regularities. 



Regularity 


r is continuous 


r is Lipschitz continuous 


solution gradient 


solution gradient 




1.8380 0.9056 


1.8650 0.9069 




1.9523 0.9513 


1.8487 0.8953 




1.7500 0.9558 


1.7570 0.9193 
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method case by case. When the solution is continuous, the loss of regularity 
has no impact on the order of convergence regardless whether the interface is 
continuous or Lipschitz continuous. When the solution is continuous, the WG 
method achieves a very good order of accuracy for continuous interface, while 
such order is slightly reduced when the interface is Lipschitz continuous. However, 
the overall order of convergence for the WG method for the latter case is still very 
close to two. Thus, it is fair to say that the WG method attains a second order 
of accuracy in L^o norm for continuous solutions. On the other hand, the Loo 
error of the WG method is always around the order of 1.75 when the solution is 
H'^ continuous, for both and Lipschitz continuous interfaces. In other words, 
even though the geometrical singularity does not affect the convergence of the WG 
method, the function singularity may do. Nevertheless, we note that a 1.75th order 
for solution with low regularities is still some of the best for this class of challenging 
problems, to our best knowledge. Moreover, although not reported in detail in this 
paper, the WG method achieves the second order of accuracy in L2 norm for the 
solution. Therefore, the present study demonstrates a good robustness of the WG 
method for solving elliptic interface problems. Readers are also encouraged to draw 
their own conclusions from the results reported in this paper. 

5. Conclusion. The present paper presents some of the best results for solving 
two-dimensional (2D) elliptic partial differential equations (PDE) with low solution 
regularity resulted from modeling nonsmooth interfaces. The weak Galerkin finite ele- 
ment method (WG-FEM) of Wang and Ye [34 is introduced for this class of problems. 
The WG-FEM is a new approach that employs discontinuous functions in the finite 
element procedure to gain flexibility in enforcing boundary and interface conditions. 
Such a strategy is akin to that used by the discontinuous Galerkin (DG) methods. 
However, by enforcing only weak continuity of variables though well defined discrete 
differential operators, the WG-FEM is able to avoid pending parameters commonly 
occurred in the DG-FEMs. Like traditional FEMs, the WG-FEM employs body-fitted 
element meshes to represent boundaries and material interfaces. Such an approach is 
convenient for handling nonsmooth interfaces or geometric singularities. The conver- 
gence analysis of the present WG-FEM for elliptic interface problems is given. The 
validity of the WG-FEM is further tested over a large number of benchmark numerical 
tests with various complex interface geometries and nonsmooth interfaces. The de- 
signed second order accuracy is confirmed in our numerical experiments for problems 
with smooth interfaces and sufficient solution regularities. 

Nonsmooth interfaces, such as sharp edges, cusps, and tips are ubiquitous both 
in nature and in engineering devices and structures. When nonsmooth interfaces 
are associated with elliptic PDEs, they result in difficulties in theoretical analysis 
and in the design of numerical algorithms. Unfortunately, to tackle the real- world 
problems, dealing with nonsmooth interfaces is non-avoidable in scientific computing. 
What aggravates the difficulty in error analysis and in the numerical methods is the 
low regularity of solution resulted from nonsmooth interfaces, which is a well-known 
natural phenomena called tip-geometry effects in many fields. When the interface is 
Lipschitz continuous and solution is or H'^ continuous, the best known result in 
the literature is of first order convergence in solution and 0.7 order of convergence 
in the gradient [21]. It is demonstrated that the present WG-FEM achieves at least 
an order of 1.75 convergence in the solution and an order of 1 convergence in the 
gradient. A numerical study for high order of WG methods should be conducted for 
a better understanding of the method. 
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